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ABSTRACT 

We study a model of of LS I +61°303 in which its radio to TeV emission is due to inter- 
action of a relativistic wind from a young pulsar with the wind from its companion Be star. 
The detailed structure of the stellar wind plays a critical role in explaining the properties of 
the system. We assume the fast polar wind is clumpy, which is typical for radiatively-driven 
winds. The dumpiness, and some plasma instabilities, cause the two winds to mix. The rela- 
tivistic electrons from the pulsar wind are retained in the moving clumps by inhomogeneities 
of the magnetic field, which explains the X-ray variability observed on time scales much 
shorter than the orbital period. We calculate detailed inhomogeneous spectral models repro- 
ducing the average broad-band spectrum from radio to TeVs. Given the uncertainties on the 
magnetic field within the wind and the form of the distribution of relativistic electrons, the 
X-ray spectrum could be dominated by either Compton or synchrotron emission. The recent 
Fermi observations constrain the high-energy cut-off in the electron distribution to be at the 
Lorentz factor of 2 x lO'* or ~ 10** in the former and latter model, respectively. We provide 
formulae comparing the losses of the relativistic electrons due to Compton, synchrotron and 
Coulomb processes vs. the distance from the Be star. We calculate the optical depth of the 
wind to free-free absorption, showing that it will suppress most of the radio emission within 
the orbit, including the pulsed signal of the rotating neutron star. We point out the importance 
of Compton and Coulomb heating of the stellar wind within and around the y-ray emitting 
region. Then, we find the most likely mechanism explaining the orbital modulation at TeV 
energies is anisotropy of emission, with relativistic electrons accelerated along the surface of 
equal ram pressure of the two winds. Pair absorption of the TeV emission suppresses one of 
the two maxima expected in an orbit. 

Key words: gamma-rays; theory - radiation processes; non-thermal - stars; individual; 
LS I +61°303 - X-rays; binaries - X-rays; individual; LS I +6r303. 



1 INTRODUCTION 

The Be star binary LS I +61°303 is one of a few currently known 
y-ray-loud (i.e., with the radiative power dominated by the y-ray 
band) X-ray binaries. The spectrum of high-energy emiss ion from 
the system extends up to TeV energies ( lAlbertetalJl2006l) and the 
peak of its eF(e) spectrum, where F is the energy flux and e is 
photon energy, is at ~ 100 MeV. 

The origin of the high-energy activity of the source is a sub- 
ject of a dispute. All binaries known to be accretion-powered have 
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either no detectable emission at e > 1 MeV or the y-ray luminosity 
much lower than the X-ray one. The latter is the case for Cyg X- 
1, where weak y-ray emission was de tected in the soft state up to 
6 < 10 MeV jMcConnell et alj|2002l) . and a weak TeV transien t 
was once detected by the MAGIC detector ^Albert et al] l2007h . 
Also, the y-ray emis sion from Cyg X-3 detected in the 0.1-100 
GeV range by Fermi l lAbdo et alj|2009bh and AGILE fravani et al] 
l2009l) is much weaker than the X-ray emission of LS I -l-61°303. 

It is possible that a y-ray loudness of an X-ray binary is related 
to its special orientation to the ob server (b y analogy with the y-ray 
loudness of active galactic nuclei, [Mirabel & Rodriguez 1999). On 
the other hand, it is possible that the binaries detected in the TeV 
y-ray band are fundamentally different from the accretion-powered 
X-ray binaries. In fact, one of the y-ray loud binaries found so far, 
PSR B 1259-63, is known to be powered by the rotation energy 
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of a young pulsar rather than by accretion ( |johnstonet"aLlll992h . 
In two other 7-ray-loud binaries, LS 5039 and LS I +61°303, the 
pulsed emission from the pulsar has not been detected, so there 
is no direct evidence for the pulsar powering the activity of these 
sources. However, the binary orbits of these two sources are much 
more compact than that of PSR B 1259-63. Then, the pulsed radio 
emission would be free-free absorbed in wind from companion star. 

If the activity of LS I +61° 303 is indeed powered by a 
young pulsar, the radio-to-y-ray emission is generated in the 
course of collision of the relativistic pulsar wind with the wind 
from the companion. Interaction of these winds leads to for- 
mation of a scaled-down analogue of the pulsar wind neb- 
ulae (PWN) in which the energy of the pulsar wind is re- 
leased at an astronomical unit, rather than on the parsec , 
distance scale 1 Maraschi & Trevei 198ll; Tavani & Aron jl 19971: 



Harrison et aT]|2 000':'s ierpows ka & Bednareldl2005l : lDubusll2006bl : 
Neronov & Che rnyakova 2007). 

In this paper, we explore the structure of the compact PWN of 
LS I +61°303 in the framework of a model with clumpy stellar wind 
from its companion Be star and assuming mixing of the pulsar and 
stellar winds. We calculate the effect of different physical processes 
that determine energy losses of high-energy particles as a function 
of the distance from the Be star. This allows us to present spectral 
models reproducing the average observed spectra from radio to 7- 
rays. The two main radiative processes are inverse Compton (IC) 
scattering of stellar photons and synchrotron emission. Coulomb 
energy loss is an important process for electrons, especially those 
with energy in the MeV range. Radio emission from within the orbit 
of the system is strongly free-free absorbed by the stellar wind. The 
relativistic electrons are carried away by the stellar wind. Overall, 
our model explains several important aspects of the behaviour of 
radio, X-ray and y-ray emission of the system. 



2 PROPERTIES OF LS I H-61°303 AND ITS WIND 
2.1 The system parameters 

The binary period o f the system, fr om radio observations, is P = 
26.4960 ± 0.0028 d l lGregorvll2003). Other binary parameters bear 
larger uncertainties jHutchings &Crampton'' 19811: ICasares et al.l 
I2OO5I hereafter C05. lGrundst'rom et al. 2007, hereafter G07). Here, 
we adopt the values given by GOtB the eccentricity of e = 0.55 ± 
0.05, the phase (0 < < 1) of the periastron of = 0.301 ± 0.01 1, 
and the mass of the Be star as M, = 12.5 ± 2.5M0. (For his- 
torical reasons, the pha se of = i s defined as corresponding 
to HJD 2,443,366.775, ICregorvllioQ^ .) In the framework of our 
model, the compact object is a neutron star, which mass we assume 
to be M2 = 1.5Mq. Then, the Kepler law yields the semi-major 
axis of a 6.3 x 10'^ cm. G07 also give the radius and the effective 
temperature of the Be star, of the BO V spectral classification, as 
R, - 4.7 X 10" cm and T, ^ 3.0 x 10* K, respectively, yielding 
the luminosity of L, = 4nRlcrT* ^ 1.3 x 10'^ erg s"', where cr is 
the Stefan-Boltzmann constant. We note that T, is relatively uncer- 
tain and some othe r studies indicate a lower value, e.g., C05 and 
IWaters et alj ( ll988l) give 2.8 x lO"* K and 2.0 x 10" K, respectively. 
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Figure 1. A schematic representation of the elliptical orbit of 
LS I +61°303, also showing interaction of the pulsai' wind with the clumpy 
Be-star wind (filled circles). The relativistic pulsar wind hits nearest clumps 
of the stellar wind. Instead of a smooth bow-shaped contact surface (shown 
by the dashed curve), the interaction region (shown by the dark circles) is 
irregularly shaped. The small ticks on the coordinate axes correspond to the 
radius of the Be star (shown by the yellow circle). The red dots along the or- 
bit mark the phase, </>, spaced with d(p = 0.05. The arrow show the direction 
of the observer, 0i„f = 0.324, <Djnf = 33°, the light shaded contour gives its 
uncertainty of ±9° , and the crosses show the superior and inferior conjunc- 
tions. The large grey circle shows the disc of the Be star. The filled red, blue 
and green squares show the points for which the orbital dependence of the 
optical depth for yy pair production is calculated, see Section|5] 



As we assume M2 = 1.5M0, the mass function of the system, 
0.011 ± 0.003Mo (C05; G07), imphes the inclination of i = 59°, 
satis fying the constraint of < 6 0° of C05, and in agreement 
with lHutchings & Cramptonl ( Il98ll) (who actually favoured higher 
values of /). Then, the separation at the periastron and apastron 
is ^ 2.8 x 10'^ cm and 9.7 x 10'^ cm, respectively. The lat- 
ter approximately equals the periastron distance of PSR B1259- 
63, and it is about 15 times smaller than the apastron distance 
of that system (assum ing the mass of its Be star of IOMq, e.g., 
iTavani & Arons Jl997!y, The distance to LS I +61°30 3 is = 
2.0 ± 0. 1 kpc I Dhawan, Mioduszewski & RupenI 20061 her eafter 
DMR06; see also lFrail & Hiellming||l99lUSteele et al.lll998l) . 

We also use the true anomaly of the system, <I> (0 < C) ^ 360°), 
which is the angle between the direction from the Be star to the 
periastron and that to the compact object. The inclination of the 
semi-major axis with respect to the line of nodes determined as = 
57 ± 9° (G07) yields the orientation of the observer, and the angular 
phases of the superior and inferior conjunctions of Oj^p = 270° - 
0) = 213° and Ojnf = 90° - a) = 33° (at the best fit), respectively. 

= 0.324. The 



^ After our calculations were completed, lAragona et al have ob- 

tained the binary parameters improved, and slightly different with respect 
to those of G07, namely e = 0.54 ± 0.03, the phase of the periastron of 
= 0.275, a) = 4l±6°, and the mass function of 0.012 ± 0.002Mo. These 
differences have negligible effect on the results presented here. 



The corresponding phases are 0sup = 0.035 and (p, 
components of LS I -l-61°303 are shown in Fig.[T| 

The mass function of the system together with the estimate 
of the Be-star mass do not allow us to determine the mass of the 
compact object, which can still be either a neutron star or a black 
hole. This is because the inclination of the orbit is only poorly con- 
strained, 10° < I < 60° (C05). The compact object would be a 
neutron star if i > 25°, and a black hole otherwise (C05). In the 
latter case, the origin of the LS I -1-61° 303 activi ty has to be accre- 
tion. T his t ype of model was first intro duced bv lTavlor & Gregory! 
l ll984h , see lBosch-Ramon et al. 1( 1200^) for a recent reference. In the 
former case, it might be accretion as well. However, the source ac- 
tivity may also be due to interactions of a young rotation-powered 



Pulsar wind nebula model ofLS I +61 °303 3 



pu lsar with the wind from the companion Be star, as first proposed 
bv lMaraschi & Trevesl l fl98ll) . 

The nature of the source has been a subject of an 
ongoing debate. Main arguments against the accretion sce- 
nario are discussed below. Its broad-band spect rum (e.g., 
IChernvakova, Neronov & Waited [20o3 : lOubusI l2006bh does not 
look like that of any black-hole accreting source at a low Ed- 
dington ratio (see, e.g., spec tra in iMcClintock & RemillartJllOOd : 
IZdziarski & Gierlinskil2004h . For example, there is no high-energy 
cut-off in the spectrum at ~ 100 keV, typical to accreting sources at 
low luminosities. Thus, attributing the spectrum of LS I -l-61°303 to 
accretion would require the existence of a yet unknown parameter 
of accretion, causing black-hole binaries accreting at the same ac- 
cretion rate to look completely different in TeV binaries and in all 
other cases. On the other hand, there is the striking similarity be- 
tween the broad band spectra of all three massive binaries detected 
as per sistent TeV so urces, PSR B 1259-63, LS I +61°303, and LS 
5039 ( lDubusll2006bl) . with the first of the objects seen to contain a 
non-accreting pulsM\ 

We note that lMassi & Kaufman Bernadol ( |2009|) have argued 
for a similarity of LS I +61° 303 to accreting black-hole binaries 
based on the value of the X-ray photon spectral index, which is 
F ~ 1.5 in the hard state of black-hole binaries as well as it 
is observed in LS I -h61°303. However, a clear difference is the 
lack of a high-energy cut-off at hard X-rays and its y-ray loud- 
ness i n the latter, as discussed above. iMassi & Kaufman Bernaddl 
( |2009|) also argued for the presence of a soft, steep power-law X- 
ray state, common in black-hole binaries, in LS I +61° 303 based 
on F = 3.6* ,^ (the Icr confidence) f ound in some INTEGRAL 
observations ^ Chemvakova et al] l2006h . However, the 2cr confi- 
dence range of those data give F between 0.7 and 10, as well as 
no strong blackbody characteristic to that state is present. Contra- 
dicting their proposal, Fermi fin ds a high y-ray flux in all states 
of the object l lAbdo et alj[2009j) . which implies that an extrapo- 
lation of the X-ray spectrum up t o 0.1 GeV has to have F < 2. 
iMassi & Kaufman Bernaddl j2009l) also notice that the index found 
in the TeV range by the MAGIC, F 2.6, is similar to the typical 
X-ray index of black-hole binaries in the soft state. However, this 
similarity is clearly accidental, and the extrapolation of that spec- 
trum down to the X-ray regime would over-predict the observed 
X-ray flux by a few orders of magnitude, again contradicting their 
supposition. 

The spectrum of LS I +61° 303 could, in principle, be domi- 
nated by the jet emission. However, the system is not viewed close 
to face-on, so the jet emission is unlikely to fully dominate the 
spectrum. Furthermore, if there is indeed a jet in the system it is 
certainly not steady, and it has been clai med to precess fast , chang- 
ing its orientation on a day time scale jMassi et"al]|2004) . Then, 
there should be periods when the jet is seen away from its axis with 
the accretion emission dominating the spectrum, but at no time the 
spectrum has been seen to resemble spectrally a standard accreting 
black-hole binary. 

The changes of the radio spatial structure of the source are 
correlated with the orbital phase (DMR06). If this is due to jet pre- 
cession, the precession period would have to be equal to the or- 
bital period, which would be a highly unusual occurrence. In all 
known cases, the disc/jet precessi on period in bi naries is much 
longer than the orbital period (e.g.. lLarwood[l998h . In the case of 
tidally-induced precession, this is also a theoretical requirement. 

If the compact object were a strongly-magnetized neutron star 
but accreting, then the system would have belonged to the class of 
accreting Be/X-ray binaries (see, e.g.. ICoell200(ll ; lNegueruelall2004l 



for reviews). However, the spectral properties of LS I +61°303 are 
very different from those of accreting X-ray pulsars in Be systems. 
Although those sources have X-ray spectra below ~10 keV often 
similar to that of LS I +61°303, they are strongly cut off in hard 
X-rays (typically with ~20 keV bremsstrahlung-li ke spectra), with 
only very weak emission above ~I00 keV (e.g.. lAppara"olll994 
iKrevkenbohm et alj2007h . 

An accreting Be/X-ray binary with the period, eccentricity and 
spect ral type similar to those of LS I +61° 303 is 4U 01 15+63A'635 
Cas jNegueruela & Okazaki|[200ll : Izi61kowskill2002h . It contains 
a B0.2 V star, P = 24.31 d, e - 0.34. If the compact object in 
LS I +61° 303 is a neutron star, then the main difference between 
the two systems can be the neutron-star spin period, = 3.61 s 
in 4U 0115+63, which is long enough to permit accretion. Also, 
its magnetic field (though strong enough not to permit equatorial 
accretion) is most likely significantly weaker, since it decays during 
accretion (Urpin & Geppert 1995; Konar & Bhattacharva 1997). 

DMR06 have strongly argued for the pulsar-wind nature of 
LS I +61° 303 based on its VLB A monitoring over one orbital 
period, showing an extended radio source with the morphology 
variable within the period and elongated along the axis connect- 
ing the stars, as expected in the case of the pulsar wind collid- 
ing with the Be stellar wind. (T his can also explai n the fast radio- 
structure variability observed bv lMassi et Z||2004 ) The morphol- 
ogy of the source varies much faster at the periastron than at the 
apastron (but no relativistic velocities have been found), suggest- 
ing the direct formation of the radio source by the p ulsar wind in- 
teracti ng with the Be stellar wind. In this context, iRomero et al.l 
( |2007|) . hereafter R07, point out that the the pulsar wind power, 
^pulsar < 10^'(SJ10'2 G)-(PJls)-\ where B, is the surface mag- 
netic field, is required to be > 10'^ erg s"' by the observed bolomet- 
ric luminosity of the system. The corresponding momentum flux is 
then significantly higher than that of the Be wind, which presents 
a problem for interpretation of the relatively coUimated structures 
observed by DMR06, as confirmed by the numerical simulations of 
R07. Indeed, we face the same problem in our calculations below, 
see Section |5] Thus, this remains to be an unresolved as yet prob- 
lem of the pulsar-wind model. On the other hand, the numerical 
simulations by R07 of their black-hole accretion model are unable 
to explain the shape of the radio structure elongated along the ro- 
tating binary axis. 

The results of t he accretion simulati ons in R07 are quite sim- 
ilar to, e.g., those in lOkazaki et al] ( l2002h . in which case the com- 
pact object is a neutron star. In fact, R07 assume both a large radius 
at which the inner disc is truncated (but they argue this does not af- 
fect their results), at which the flow structure would not be affected 
by the nature of the central compact object, and a low black-hole 
mass of 2.5Mq, also not very different from the possible masses 
of neutron stars. Yet no accreting neutron-star binary has been ob- 
served to be similar to LS I +61°303. 

R07 also argued that the lack of pulsations detected from 
LS I +61° 303 represents an argument against the pulsar model. 
However, as we show in Section |273l below, this lack is actually 
expected given the compactness of the orbit. 

A further evidence for the presence of a young pulsar is pro- 
vided by a detection of a powerful ~0.2-s flare with the peak 15- 
150 keV flux of 5 X 10"* erg cm"- s"' (which is higher by a factor 
of ~ 10^ than the avera ge source flux) and a blackbody-like spec- 
trum by the Swift/BAJ' jPe Pasguale et al.ll2008l : iBarthelmv et all 
l2008l) . Such an event is consistent with a soft gamma-r ay re- 
peater/anomalous X-ray pulsar burst l lDubus & Giebelsl2008h . This 
would indicate the presence in LS I +61°303 of a very strong mag- 
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netic field, typ ical to magne t ars. S ome of those objects do emit pul- 
sar wind, e.g.. lGavriil et alJjlOOSh . as required in the present case. 

An argument for the black-hole nature of the compact object 
(thus implying accretion as the underlying cause of the source ac- 
tivity) is provided by some measurements of the width of absorp- 
tion lines in the system. C05 and GOV have found that width to 
be vsin; 100-110 km s"'. It is most likely due to the rotation 
of the Be star. Be st ars usually ro tate at ~0.7-0.8 times the criti- 
cal rotation velocity l lPorteij|l996i) . which implies ; ^ 15°-20° in 
the case of LS I -l-61°303 though a lower rotation velocity of the 
star cannot be ruled out (COS). This inclination would then im- 
ply the compact object is a black hole. However, G07 argue that 
those measurements are affected by weak (and time-dependent) 
emission in the line wings that made the absorption cores appear 
more narrow. This then can r econcile these results with those of 
iHutchings & CramptonI Jl98ll) . who obtained v sin / 360 ± 25 km 
s"', and whose results imply i > 50° (M2 < 2M0), in full agree- 
ment with the neutron-star nature of the compact object. 

We also note that all Be binaries with know n com- 
panions contain neutron s tars (e.g., iNegueruelj |2004 
iBelczvhski & Ziolkowskil I2OIOI) . (White dwarf and He subd- 
warf companions are also predicted, but are very difficult to find.) 
The reason for that may lie in the system evolution. The rotation, 
being the defining characteristic of Be stars, may be achieved 
during a period of Roche-lobe overflow mass transfer from its 
(initially) more massive companion onto a B star. However, the 
companion loses most of its mass in the process, becoming a He 
star with the mass of only a few Mq, which, when explodes as a 
supernova, can l eave behind only a neutron star, but not a b lack 
hole (see, e.g., iGiesI 120001 : | Tauris & van den Heuvell [2OO6I and 
references therein). On the other hand, if rotation of Be stars is 
ach ieved in a different way, the popu lation synthesis calculations 
by IBelczvhski & Ziolkowskil j2010l) show that Be/black-hole 
binaries are evolutionary formed only in ~ 1/30 of cases compared 
to those with neutron stars, which can explain the observed absence 
of black-hole systems. Still, the existence of < 2 black holes in 
known Galactic Be X-ray binaries is possible. 

2.2 Structure of the stellar wind 

Be stellar outflows have generally two components. One is 
a radiation-driven and approximately isotropic, fast wind. The 
other is a slow equatorial outflow forming a thin decre- 
tion disc. The pr ocess causing the form ation of the disc 
remains unknown jPorter & RiviniusI l2003i) . The disc seems 
to be viscous and nearly Keplerian, with the radial veloc- 
ity much smaller than the rotational one. Observationally, 
Be decretion discs in binaries ha ve outer radii significantly 
smaller than those o f isolated stars l lReig. Fabregat & Coelll997l ; 
IZamanov et al.l I2OOII ). This appears to be caused by tidal trun- 
catio n (via Lindblad resonances) of the disc by the neutron 
star Artymowicz & Lubow' 'l994'; 'Negueruela & Okazakij 1200 ll ; 
lOka zaki & Negueruela 2001; Okazaki et al. 2002). Numerical sim- 
ulations specific to LS I +61° 303 showing tidal truncation of the 
decretion disc have been performed by R07. 

The size of the decretion disc can be measured via the equiv- 
alent width of the Ra line, W,i. In LS I +6r303, varies be- 
t ween -6A and -18A (GOT). Using the numerical Be disc model 
of iGrundstrom & GiesI ( |2003) and assuming ; = 60°, this yields the 
disc radius varying between 4R, to 6R,. However, COS found that 
in order to match the spectrum of the BO V star in LS I -l-61°303 to 
a template (of the BO V star t Sco), the template needed to be re- 



duced in the flux. If the resulting additional stellar flux comes from 
the decretion disc (as postulated by COS), W,\ needs to be increased 
(to account for the stellar continuum being lower) by 1.54. This (a 
relatively uncertain) correction has a minor effect on the disc ra- 
dius, increasing it to (4.5-l)R,. On the other hand, the separation 
varies in the (6. l-2I)if, range (Section |2j. Thus, the neutron star 
can directly interact with the disc only around the periastron, and 
only when the circumstellar disc reaches its largest size. On the 
other hand, the disc is not circular, and some disc matter can reach 
the neutron star around the periastron even if the average radius is 
le ss than the separation , as shown in the hydrodynamic simulations 
of lOkazaki et alj ( |2002|) . 

The density profile of the equatorial disc in LS I +6I°303 
ha s been measu r ed usi ng t he IR free-free and free- bound radiation 
bv I Waters etalj jl988l) and lMarti & Paredesl h995l) . Assuming the 
disc half-opening angle of 9o = 15° and R, = 7 x 10" cm, they 
found the disc ion density of 



(1) 



where «d.o ~ 10'^ cm"^, and q ^ 3.2. Here, D is the distance 
from the center of the Be star, which we hereafter will take as a 
characteristic distance. The disc density at the periastron separa- 
tion (D ^ 3 X 10'^ cm) is then ~ 3 x 10'" cm"^ The actual value 
of iiiifi is relatively uncertain, and it will be somewhat different at 
our adopted R,, as well as it depends on the unknown value of 9q. 
For example. IPorted ( fl996h found that {Og} in Be stars is 5°. Also, 
rii varies during a long-term cycle of the disc activity. 

Furthermore, it is highly uncertain to what mass-loss rate 
equation l[T]l corresponds. This depends on the initial outflow ve- 
locity, Vdo, which has n ot been measured. I Waters et al.l l ll988l) and 
iMartf & Paredesl l ll995h assumed 5 km s ' and 2-20 km 



spectively. However, Ported ( Il998l) found that such outflows would 
result in a fast removal of the angular momentum from Be stars, and 
their substantial spin-down, which is inconsistent with rotational 
velocity distributions for different Be luminosity classes, showing 
no significan t evolution. Thus, likely va lues of i'd,o are « 1 km s"'. 
For example. lMarlborough et al.l ( [l997l) used Vd,o - 0.3 km s ' (and 
9q = 5°) as the standard disc parameters. At this Vd.o, the equatorial 
mass-loss rate is Mi ~ 10"** Mq yr"', and the outflow velocity at 
the periastron distance, implied by equation ([T), is ~3 km s"'. 

The other outflow component is a radial, radiation-driven, 
wind. The wind velocity for such outflows follows the usual law. 



VwiD) - I'w.O + (Voo - I'w.o) 



('-I) 



(2) 



where Voo - (l-2)x lO' km s"' and Vw.o ^ are the initial and ter- 
minal, respectively, wind velocity, and /? ~ 1. The number density 
averaged over possible wind clumping (see below) at D is related to 
v„(Z)) and the wind mass-loss rate, M„, via the continuity equation. 



<«»,(£•)> ^ 



47tmpyUi£)2vw(£>) 4nmpii,D(D - R,)vo 



(3) 



where yUj ^ 4/(1 + 3^) is the mean ion molecular weight, which 
is 1.3 for the H abundance of 0.7, and the second equality is 
for/? = 1 and neglecting \'„_o- Unfortunately, the wind parameters 
have not been measu red for LS I -l-61°303. Following the general 
estimates of W aters et al.l ( Il988l) . we use fiducial values of M„ = 
10"* Mq yr"' and Voo = 10^ km s"'. The kinetic power of the wind 
isMwV^/2^3x 10" ergs-i. 

Radiation-driven winds from massive stars are observed to 
be clumpy, with the clumps filling a fraction of / ~ 0.1 of 
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the wi n d volume, see, e.g., Moffat & Robert 1 1994 ). Schild et all 



ilOOA liouret. Lanz & Hillieil bOOSi) . IPuIs et al.l ( l2006h . Formal 
tion of clumps in these winds is predicted theoretically, with hydro- 
dynamical instabilities leadin g to formation of shocks, density en- 
hancements and rarefactions dOwocki. Castor & Rvbickilll988h . In 
the case of Be stars, density irregularities on the scale of ~ 10" cm 
were found in the wind of 2S 0114-1-65 jApparao. Bisht & Singhl 
Il99lh . The wind density inside the clumps is 1 // times that of the 
smooth wind. For the parameters adopted above and at D » i?,, we 
thus have. 



(4) 
(5) 



\3 X 1012 cm/ 



\ lO-^M© yr-' 



The presence of clumps affects the interaction with the pulsar wind. 
Also, the rates of two-body processes, e.g., free-free emission and 
absorption, are enhanced by a factor of 1 // with respect to the case 
of a smooth flow. The detailed form of clumpy winds remains un- 
known; in particular, it is likely that the space between the clumps 
is not empty. In that case, the effect of clumping on a process with 
the rate proportional to the square o f a density, n, can be des cribed 
by the generalized clumping factor jOwocki & Cohenll2006h . 



(6) 



where the averages are along a line of sight within a volume of ap- 
proximately constant properties. If the density is constant within 
the clumps and the space outside is empty, this / becomes the 
usual volume filling factor We note that when individual clumps 
become optically thick, porosity of the clu mpy wind becomes 
important and reduces t he effective opacity jOwocki et al .11 20041 ; 
lOwocki & CohenlbOOd ; lOwocki et al.ll2009<) , which complication 
we neglect in Section l231 below. 

Apart from the density, the other main parameter of the wind 
is its temperature, T . In the case of a single massive star, T cor- 
responds to radiative equilibrium in the photon field of the star. 
It is then somewhat l ower than T,, and in the case of Be discs, 
ICarciofi & BiorkmarJ find T ~ 0.6T,. For the fast wind, it 

is likely to be somewhat higher. The wind temperature reaches the 
Compton temperature of 0.967'» in the fuU-ionization regime. 

However, X-ray and soft y-ray emission substantially in- 
creases the (optically-thin) wind temperature even if its luminosity 
is much lower than the stellar one, because the Compton heating 
rate (in the Thomson regime) is proportional to ^ EF{E), where 
F{E) is the energy flux of a local spectrum. This effect leads, e.g., 
to an increase of the temperature of a part of the wind close to the 
X-ra y source to values as high as ~ 10 ^-10^ K in Cyg X-1 and Cyg 



X-ra y source to values as high as ~ 10 
X-3 Iszostek & Zdziarskill2007l . l2008h 



We thus calculate the Compton temperature in LS I -1-61° 303. 
Unlike the case of accreting sources, most of the Compton heat- 
ing here is done by y-rays, not by X-rays. Then, the calculations 
need to be done employing the full Klein-Nishina cross section 
rather than using the usual Thomson-limit method (used, e.g., by 
[Sege lman. McKee & Shields 1983)- We have calculated the net 
rate of energy transfer betw een thermal elec trons and a given pho- 
ton field using the results o f iGuilberll ( fl986l) . We use the X-ray/y- 
ray broad-band spectrum as given in Section 13.31 which can be 
roughly approximated as a doubly-broken power law. 



eF{e) 



1 erg cm"- s 



1.6x 10-"(j|v) , e<lMeV; 

1 X 10-'"(Yf^)°, IMeV < 6 < lOGeV; (7) 

, / \-0.76 

4.1x10-3(^1^) , e>10GeV 



(corresponding to the high X-ray state), with the luminosity of 
Lxy - 6 X 10^^ erg s"' (assuming isotropy and d = 2 kpc). In 
addition, there is the stellar blackbody at = 3 x 10* K and 
L, = 1 X 10^^ erg s"'. We find the Compton temperature for the 
composite spectrum is T = 4 x 10^ K, i.e., more than an order 
of magnitude higher than T,. Heating by photons with energies > 1 
MeV dominates that by photons below 1 MeV, in spite of the Klein- 
Nishina decline of the Compton cross section. 

The above composite spectrum corresponds to points at the 
same distance from each of the two sources of radiation, and ne- 
glects the presumed diffuse nature of the high-energy source. The 
actual relative proportions of the high-energy spectrum and of the 
blackbody depend on the position. In a region around the y-ray 
source, the temperature will be substantially higher than the T 
above. In particular, the Compton temperature for the high-energy 
spectrum alone, equation Q, is as high as 7 x 10** K. On the 
other hand, line cooling will strongly reduce the temperature away 
from the high energy source (compare the detailed wind structure 
in Cyg X- 1 , Szoste k & Zdzim^ski 200Z), which effect increases with 
the clumping (e.g., ISzostek & Zdziarskill2008l) . 

The Compton cooling is almost entirely du e to stellar pho- 
tons, and its rate per electron can be written as (e.g.. lBegelman et"ail 
Il983h . 



C 



kT a-jL, 



(8) 



1.6 X 10" 



'( ^ )f ]( ^ \ 

Ux 105K/\l038ergs-VV3 X 10'^ cm/ 



erg s 



where k is the Boltzmann constant, tr-p is the Thomson cross sec- 
tion, and nie is the electron mass. In equilibrium, it equals the heat- 
ing rate. This can be compared with other relevant rates. Heating 
of the wind by Coulomb interactions with relativistic electrons (see 
Section [J!4l below) may play a major role here. We find that heating 
rate, Lc/N (where A' is the total number of particles in the heated 
volume, ~ 2£)'(«i)), to be comparable to the rate of equation ^ 
already when the total power supplied to the wind by the Coulomb 
process is Z-c ~ 10" erg s"', i.e, a ~ 2 x lO^^Lxy. If Lc ~ lO^"* 
erg s"', the equilibrium (but neglecting line cooling) wind temper- 
ature will be ~4 X 10^ K. On the other hand, a major cooling effect 
is due to advection of the thermal energy by the wind away from 
the central region, ~ IkTv^jD. At Vw - 10** cm and 4 x 10^ K, 
this rate is ~ 4 x 10"'^(D/3 x 10'^ cm)"' erg s"'. These effects will 
strongly depend on position via the spatial dependence of both the 
wind density and of the relevant rates. 

These details are beyond the scope of this work, intended to 
identify the main relevant physical processes. Given the consid- 
erations above, we use a fiducial value of the wind temperature 
of r = 10^ K. We note that the measurements of the equatorial 
mass loss rat e discussed above do no t depend on T (apart from the 
Gaunt factor. I Wright & Barlowll 19751) . On the other hand, free-free 
absorption strongly depends on T, see Section l231 below. 



2.3 Free-free absorption 

The free-free absorption coefficient due to ions with the atomic 
charge, Z, is given by (e.g. lRvbicki & Lightmanlll979l) . 
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Q-ff 



25/2„l/2^ 



{kT)-'''Z'n,nzV-'g, 



(9) 



where v is frequency, «z is tlie Z-ion density, is the electron den- 
sity, e is the electron charge, and g is the average Gaunt factor The 
Gaunt factor for hv <k kT and v » Vp equals 



31/2 



In 



(10) 



wher e Vp is the pla sma frequency and ^ 0.5772 is Euler's con- 
stant ( ISDitzeilll978h . Then, using Z = 1, T = 10' K and v = 5 GHz 
in equation ^W^, and averaging over Z in equation ([9}, we find 



(11) 



where T and y are in units of K and Hz, respectively, and = 
2/(1 -I- X) 1.2 is the mean electron molecular weight. 

We first calculate the optical depth of the equatorial disc, 
which density is given by equation l[TJ. We consider the optical 
depth perpendicular to the disc plane. The disc thickness equals 
2Z)tan6o. Thus, 



Td.ff 



(12) 



^ U0i3cm-3/ llOSK/ llGHz/ UxlO'^cm/ ' 

Thus, the entire equatorial disc is completely optically thick to ra- 
dio emission. 

Then, we calculate the radial optical depth of the fast wind 
from infinity down to a given value of D. In a clumpy medium, 
Tff, is an integral over the square of the density within the clumps 
times /, i.e., n^f (or equivalently, an integral over (/J;)"/"'). Using 
equation (|4j, we obtain. 



^"•f " ^ llO-«Moyr-i) ( 10« cm s-') (ot) 

^ (iGHi) (t05k) (3 X 1012 cm) ' ''^^^ 



which corresponds to Tw.ff = 1 at 



5 X 10' 



2/3 



10-*Mo yr-' 

2/3 , J v-l/2 



\10«cms-i/ \0.1 



UGHz/ VIO^K/ 



-1/3 



(14) 



Our results imply that (a) the neutron star moves entirely in 
the optically thick region, explaining the absence of observed radio 
pulsations and (b) that the radio emission from the system has to be 
produced at the distances D > D^ g, lying much outside the binary 
orbit. The optical thickness of the cent ral region of LS I +61°3() 3 
has been pointed out before, e.g., by iTavlor & Gregorvl ( Il982h . 
iMassi et al.l ( 1 19931) and lDubuj f2006bf) . An additional process ab- 
sorbing radio photons in the central re gions is synch rotron absorp- 
tion, which is likely to be substantial ( lLeahvll2004h . On the other 
hand, the ~ 1 MeV electrons, responsible for the synchrotron emis- 
sion from the region of interaction of the pulsar and stellar winds, 
see Section l?!4l below, can escape beyond the distance £)ff, so that 
the radio emission from the extended wind interaction region can 
be much less suppressed than the direct pulsed emission from the 
neutron star 



3 WIND INTERACTION 

3.1 Homogeneous model of the wind interaction 

In the simplest version, the model of interaction of the pulsar wind 
with the wind from the companion massive star assumes that an 
isotropic relativistic outflow from the pulsar hits a smooth (but 
anisotropic, in the case of a Be star) outflow from the star. Geom- 
etry of the interaction surface is bow-shaped, as determined by the 
pressure balance between the pulsar and stellar winds. The physi- 
cal parameters of pulsar winds remain poorly known. In a simple 
model, an isotropic e* wind with a fixed bulk Lorentz fact or is as- 
sumed . Such a model was developed in details b^ Tavani & Aronj 
( Il997h . and applied to LS I -^61=303 bv lDubuj ( l2006bl) . 

In this model, the pulsar and stellar winds do not mix, with the 
boundary between the shocked and un-shocked stellar wind being 
spatially separated from the boundary between the shocked and un- 
shocked pulsar wind. This allows the shocked pulsar wind to escape 
with a speed of ~ c/3 (which is much higher than the velocity of 
the shocked stellar wind). The corresponding time scale is. 



:^.30of — ^^)s, 

c/3 V3 X lO'^cm/ 



(15) 



where is the characteristic smooth wind size (taken as the pe- 
riastron separation). However, the power of the pulsar wind itself 
is constant. This means that any variability is determined by the 
details of the pulsar wind interaction with the stellar wind rather 
than by the above time scale. The system properties change on the 
orbital time scale. 



^ ~ 300( ^ \{ X 

i'o,-b \3x 10'2cm/V10'cms-i/ 

where Voa is the pulsar orbital speed. 



ks. 



(16) 



3.2 Inhomogeneity of the stellar wind and short time-scale 
variability 

A long XMM observation of LS I -l-6r303 in 2005 has revealed 
variability of the system at the time scale. At, of several ks 
l lSidoh et al.ll2006l) . which is much shorter than the orbital period. 
Variability at a similar time scale i s also observed in the radio band 
jPeracaula. Martf & Paredeslll997l) . This implies the presence of a 
variability mechanism much faster than overall changes of the sys- 
tem properties with the orbital motion, equation 1 II6I 1. This is pos- 
sible if the characteristic size scale of variability of the stellar wind 
is much smaller than the binary separation. This implies that the 
stellar wind is clumpy rather than smooth, as we have argued on 
independent groun ds in Section|2j2] Furthermore, recent results of 
ISmith et al. 1 120091) have also shown rare events in the form of flares 
with the doubling time scale as short as > 2 s. 

If the wind is clumpy, then the observed variability time scale 
of ~ 10 ks is likely to correspond to the time of the passage of a 
clump through an interaction region. Using v„ ~ 10^ cm s"' of the 
polar wind, the smaller of the required clump size, and the size 
of the interaction region is ~ i'„A? ~ 10'2(v„/10* cm s"')(Af/10ks) 
cm. However, the characteristic size of a clump can be smaller if 
the observed variability results from averaging over passages of a 
number of the clumps. The > 2-s time scale may correspond to 
very rare events when the pulsar wind manages to irradiate a few 
clumps with a high magne tic field, as expected close to the stellar 
surface. ISmith et aU ( |2009|) also discuss the possibility of this time 
scale corresponding to the stand-off distance between the pulsar 
and stellar winds. 
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Figure 2. (a) The radial profiles of the stellar-wind density and velocity 
(each normalized to unity at the maximum) assumed in the spectral calcu- 
lations shown in Fig. [3] The resulting normalized flux radial profiles [per 
ln(Z))] are shown in (b), (c) for the cases computed in Fig. [31 a) and (b), re- 
spectively. The vertical dashed line marks Do. the assumed lower boundary 
of the energy injection region. The black long-dashed curve in (b) shows the 
fraction of energy of 10-MeV electrons going into Coulomb losses. There is 
no 1-GHz curve in (c) because that model yields virtually no radio emission. 



Furthermore, the clumpy stellar wind and the pulsar wind are 
very likely to mix, e.g., due to plasma instabilities. Then, magnetic 
field can retain the high-energy particles within the clumps, as we 
show in Section [34l below. Then, their escape velocity slows down 
to the stellar wind velocity. The corresponding escape time is. 



£)/v„ ~ 3 X 10** 



D 



,3 X 10'2( 



i)(lO'* cms-') 



(17) 



It is not clear whether the multi-TeV electrons responsible for 
the observed TeV y-ray emission are initially present in the un- 
shocked pulsar wind, or they are produced via shock acceleration 
in the wind collision. Another mechanis m of injection of relativistic 
electr ons can be proton-proton collisions jNeronov & Chemvakoval 
l2007h . Here, we thus assume an injection of relativistic electrons 
and study its consequences in the presence of wind mixing. 



3.3 Spectral models 



from superposition of con- 
this, we provide a general 



Below we calculate spectra resulting 
tributions from difl^erent radii. Before 

characterization of th e source in terms of the photon and magnetic 
comp actnesses (e.g.. IVurm & Poutanen|[2003 : [ 
l2009l) and the Thomson optical depths. 



Malzac & BelmonJ 



tn = 



IQ- 



{B-/Sn)o-jD 



\(fi^ erg S-' 

^Xy 
1036 gj-g g-l 

B 



These compactnesses are, 
D 



10" 



V3 X 1012 cm/ 
V3 X 1012cm/ 



\1G/ \3 X lO'^cm/ 



(18) 



(19) 



(20) 



where B is the magnetic field strength. Thus, all the compactnesses 
are ^ 1, and Ib ^Xy <K ^ 1. In the Thomson regime, single 
Compton up-scattering of the stellar photons is usually the main 
radiative process. However, electrons with the Lorentz factor of y > 
10' are in the Klein-Nishina regime, and then synchrotron losses 
may dominate. The magnetic field is typically below equipartition 
with the radiation energy density. 

The equatorial Thomson optical depth in the disc component 
of the wind, equation lUll, from D to infinity (neglecting the disc 
truncation) is given by. 



Td ^ 0.02 



\10i3cm-3/\3 X 1012cm/ 



(21) 



The radial Thomson optical depth in the radial component of the 
wind, equation ([3]l, from D to infinity (neglecting its truncation) is 
even lower. 



6x 10''( ,n8'°° i ) 
\ 10'* cm s-i / 



D 



10-8Moyr-i;\3 X IO12 cm 



.(22) 



Thus, the Thomson optical depths of interest are « 1. 

We have constructed a one-dimensional spectral model of the 
source, taking into account the synchrotron and bremsstrahlung 
emission, free-free absorption by the stellar wind, and Compton, 
Coulomb, and photon-photon e* pair production processes. We 
consider a region of the wind in the equator moving with the wind 
velocity (i.e., assuming that the relativistic electrons are trapped in 
the wind, see Section [T4l below). Then, the time evolution of the 
electron spectrum corresponds also to the spatial movement. As our 
one-dimensional model is taken to represent the entire source, the 
considered volume corresponds to a conical region centred on the 
star. We thus use the standard kinetic equation for the evolution of 
a relativistic electron distrib ution in the presence of radiative losses 
l lBlumenthal & Gould|[l970l) , in which we replace the time deriva- 
tive, (9/3/, by the spatial one, v{D)dldD. The escape of photons 
is also treated one-dimensionally, outward along the radial coordi- 
nate. 

We follow the evolution of electron spectrum in the distance 
range of Z)o = 2 x 10i2cm < £) < lOi"* cm. The power in the 
injection per unit radius in the considered conical volume is as- 
sumed to be oc The electron injection spectrum is assumed to 
be an e-folded power law, AN^IAE oc (E + m^c^Y^^ exp(E/Ec) at 
£ > 0, where E is the electron energy, and E^ is its cut-off value, 
E^ = jcm^c^. The pow er-law injection spectrum may originate 
in the shock region (e.g.,lBogovalov & Aharonian''2000!), resulting 
from the wind collision (e.g.. ,Khangulvan et al..2007i) . 

The radial profiles of the density and velocity of the wind as- 
sumed for the spectral calculations are shown in Fig. [3 a). The ac- 
tual wind contains the disc and radial components, see Section |T2l 
which would be difficult to take into account in one-dimensional 
calculations. Therefore, close to the star, at £) < 6R,, we consider 
only the dense circumstellar disc, with the density following equa- 



tion lO 



assuming rii^o 



X 1012 cm"^ and q = 3.2. The initial 



Then, the disc is truncated at D = 6R, and outside it the wind 
asymptotically reaches the terminal velocity of the radial wind of 
Voe = 10^ cm s^i. Thus, the density decreases at large distances as 
n(D) oc D^^, see equation l[3]l. For simplicity, we neglect in these 
models the effect of the wind dumpiness on the radiative rates. 

The ma gnetic fields of Be stars remain mostly unknown 
l lNeinedl2007h . though dipol e magnetic fiel ds of several hundred G 
were claimed in some cases jStepiehil2003i and references therein). 
On the other hand, the upper limits on the Be-star fields found by 



8 Zdziarski, Neronov and Chernyakova 




-5 -4 -3 -2 -1 1 2 3 4 5 6 7 8 9 10 11 12 13 

log e[eV] 



Figure 3. Model s pectra comp ared to the dat a, which are the same as those in IChemvakova et alj fcOOd) except for the data from the VERITAS 
jAcciari et al.ll2008h and Fermi jAbdo et alj|200 9a) instruments. The dot-dashed, dotted and dashed curves show the spectral components from the syn- 
chrotron, bremsstrahlung and IC processes, respectively, and the thick solid curve gives the sum. (a) The model in which the X-rays are dominated by the IC 
process, (b) The model in which the X-rays are dominated by the synchrotron process. See Section l33] for the parameters. 



lHubrigetal] ( l2007h are at the level of tens of Gauss. Then, a dipole 
spatial dependence would imply the field at Do of < 1-10 G. At 
D > Do, we assume B = Bq(Do/D), and adopt Bo = 2.3 G. 

Fig. |3] shows the results of calculations for two choices of 
the electron injection. Figs.[3la), (b) show results for the cases of 
Te = 1.7, £, = 10 GeV (y, ^ 2 x 10"), and T, = I, E, = 50 TeV 
(7c - 10**), respectively. In the two cases, the dominant contribution 
to the 1 keV-10 GeV spectrum is the IC and synchrotron emission, 
respectively. Both m odels fit well the Fermi 0.1-100 GeV spectrum 
jAbdo et al.ll2009al) . in particular, its relatively sharp high-energy 
cut-off. The bolometric luminosities are 5.4 x 10^' erg s"' and 
5.1 X 10^^^ erg s"', respectively [isotropic at d = 2 kpc, and in ap- 
proximate agreement with the approximation of equation l|7j], and 
the contribution of the non-thermal bremsstrahlung is minor in both 
cases. 

In the case shown in Fig. ^a), the high-energy cut-off in the 
electron distribution at yields a spectral cut-off in the Compton 
spectrum reproducing the Fermi data at e > 3kT,yl 3 GeV. An 
additional minor effect increasing the sharpness of the cut-off is 
the onset of the Klein-Nishina regime of Compton scattering, at 
e > (m^c-f- 1 {?>kT ,) ^ 30 GeV or so. We see that the model with the 
dominant IC emission fails to explain the TeV emission, which then 
needs to be accounted for by a separate spectral component. On the 
other hand, preliminary results of observations of LS I -1-61° 303 
by VERITAS contemporaneous with those by Fermi have yielded 



no detection, with the upper limits on the 1 TeV flux a factor of 
~ 2 below the fluxes from MAGIC and VERITAS at other epochs. 
Thus, the TeV emission during the Fermi observations was at most 
weaker than that shown in Fig. [3] 

The steady-state electron distribution, dN^/dE, is determined 
by the form of the injection and the energy loss and escape pro- 
cesses. Then, the break in the bremsstrahlung spectrum around 
e ~ 30 MeV is due the transition from the dominance of Coulomb 
losses at E\, < 30 MeV (yb - 60) to the dominance of Compton 
losses at higher electron energies. In the former regime, roughly 
dN^ldE cc dN^/dE, while in the latter, dN^/dE is steeper by unity, 
cx y"^ '. This break is also reflected in the Compton spectrum, at 
e ~ 3kT,jl ^ 30 keV. 

We see that the synchrotron component of this model accounts 
quite well for the radio flux. The peak of this spectrum corresponds 
to characteristic energy of synchrotron photons emitted by an elec- 
tron with 7c, e - Q/4n)yl(heBo/2m^c) ~ 10 eV, where h is the 
Planck constant, e is the electron charge, and 1/2 accounts for the 
average sine of the pitch angle. The part of the steady-state elec- 
tron distribution emitting the synchrotron emission is an e-folded 
power law, and then can be calculated analytically (see eq. 9 in 
IZdziarski etai]|2003h . The break at e > 0.01 eV (y > 2 x 10'- Hz) 
corresponds to the onset of free-free absorption by the stellar wind, 
see Section l23l 

In the case shown in Fig. (S^b), the spectral component dom- 
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inant at e < 20 GeV is synchrotron emission. Its dominance over 
Compton scattering is due to the electron injection being strongly 
concentrated at ^ 10*, which is very deep in the Klein-Nishina 
regime (which corresponds to 7 > 10^). This strongly reduces the 
Compton losses and allows the synchrotron losses to dominate at 
electron energies responsible for most of the synchrotron emis- 
sion. The peak of it is close to the characteristic energy of e ^ 
(3/4n)jl(heBo/2m^c) ~ 200 MeV (see above). The steady-state 
electron distribution is approximately oc y"^ at both the regimes 
dominated by synchrotron losses and by Thomson losses (7 < 10^), 
but at different normalization in each regime, with the regime dom- 
inated by Klein-Nishina losses in-between. 

Then, the Compton component dominates at photon energies 
at which synchrotron emission is by electrons deep in the exponen- 
tial tail of their distribution, > lOyc, e > 20 GeV. The Compton 
component breaks at the transition from the Thomson to Klein- 
Nishina scattering, 6 ~ 100 GeV, corresponding to 7 ~ 10'. This is 
also the onset of substantial pair absorption on the stellar photons. 
On the other hand, the present model cannot account for the radio 
emission, which needs to be due to a different physical component. 
Given that the previous model accounts for the radio emission but 
not for the TeV one, the actual physical conditions in the source 
may correspond to a superposition of the two models considered 
here. 

An important issue is the origin of the very hard injection 
and the high maximum Lorentz factor, =^ 10** in this model. 
Very hard, almost mono-energetic acceleration rates have been pre- 
dicted theoreticall y and invoked to explai n some astrophysic al phe- 
nomena, see, e.g., iLesch & Reici] ( ll992l) . |Protheroel ( 12004) . How- 
ever, a problem is presen t ed by the high value of y^.. As noted by 
iGuilbert. Fabian. & ReesI ( Il983h . the acceleration rate is generally 
oc B, the synchrotron energy loss is oc S", and the synchrotron pho- 
ton energy is oc By-. Combining these dependencies lead to the 
maximum possible characteristic synchrotron photon energy from 
accelerated electrons being independent of B, 

£max,s = -| — - 20^ MeV, (23) 

32 Q-f 

where ^ < 1 is an acceleratio n efficiency and at is the fine-structure 
constant (cf., e.g., eq. 8 in IZdziarski. Malzac & Bednarekl [2009I . 
where, however, no average over the pitch angle was included). On 
the other hand, the characteristic synchrotron energy in the spec- 
trum in Fig.[3tb) is =^ 200 MeV, 10 times higher even if we take the 
maximum efficiency of ^ = 1. Thus, a dilferent mechanism than 
shock acceleration is needed. This can be due to direct injection of 
electrons from cold relativistic pulsar wind with the bulk Lorentz 
factor equal to 7c. Alternatively, the elect rons could be produced in 
interactions of high-energy proto ns (e.g.. lNeronov & Chernvakoval 
l2007l : lNeronov & Ribordvil2009l) . 

Figs.|2lb-c) show the normalized radial flux profiles at photon 
energies in the radio, X-ray and high-energy 7-ray regimes, corre- 
sponding to the models shown in Figs.[3ta-b), respectively. We plot 
here the flux per logarithm of the distance. We see that whereas the 
X-ray and 7-ray fluxes are dominated by the contributions from 
around Do, the radio flux in the model of Fig. [S^a) peaks above 
10'^ cm. We also show the flux profiles in a different representation 
in Fig.|4l in which the actual fluxes are plotted in such units that 
they are proportional to the surface brightness. 

We have also quantified the role of Coulomb losses of rela- 
tivistic electrons for the model of Fig.[3la). First, we have calcu- 
lated the fractional energy loss going into Coulomb losses for 10- 
MeV electrons, shown by the black long-dashed curve in Fig.|2{b). 



It is > 0.5 at £) < 6 X 10'^ cm. Then, the corresponding profile of 1 
keV emission would be very strongly centrally peaked if Coulomb 
losses were not included. We have also calculated the spectrum as- 
suming no Coulomb losses (not shown here), in which the effect 
of Coulomb losses is seen for I eV < e < 10 MeV. That spectrum 
is approximately flat [in eF(eJ\ already at e > 10 eV. The power 
going into Coulomb heating of the wind in this model is very sub- 
stantial, 1.5 X 10^^' erg s"'. On the other hand, no such effects 
are present in the case of Fig.[3tb) because both the shown model 
spectrum is produced by electrons with 7 > 2 x 10' as well as the 
contribution of low-energy electrons to heating is negligible due to 
the hard injection. 

At e > 10 GeV, the shape of the photon spectrum starts to 
be modified by photon-photon e- pair production. The soft pho- 
tons responsible for the pair opacity originate from both the Be 
star and its circumstellar disc. However, such discs are transpar- 
ent in the optical range (since no partial eclipses of Be stars by 
discs are observed), and thus their emission has to be much weaker 
than the blackbody spectrum at the disc temperature (^ 0.6T, 
following ICarciofi & BiorkmanI |2006|) . Indeed, detailed calcula- 
tions of Be disc emission show a substantial excess of emission 
above that of the star only at e < 0.5 eV (^ 2^m), see, e.g., 
ISigut. McGill & Jones! ( l2009h . This also roughly agrees with our 
calculations in Section 12.31 which show that the disc at Dq be- 
comes optically thin to free-free absorption at y 3 x 10" Hz 
(^ O.I eV), i.e., the emission is blackbody only up to this en- 
ergy, and becomes quickly much weaker at higher energies. Thus, 
we approximated the disc emission as Rayleigh- Jeans at < 0.1 
eV and as optically-thin bremsstrahlung (neglecting the free-bound 
contribution to emission) at > 0.1 eV, with the disc truncated at 
6R,. Though photons at 1 TeV have the threshold for the e* pair 
production at ^ 0.3 eV, the stellar blackbody spectrum increases 
fast with energy and the pair opacity is dominated by photons at 
higher energies. Thus, absorption by disc emission has a relatively 
minor effect on the shape of the observed TeV emission, except 
ab ove a few TeV. We notice that our approach differs from that 
of lOrellana&Romerol l l2007h , who assumed the disc emits only 
slightly diluted blackbody at the disc temperature. This disagrees 
with the Be disc physics, and would strongly overestimate the 7- 
ray absorption by its emission. Then, R07 have calculated the de- 
tailed disc emission in their model, and took it into account for pair 
absorption. However, they do not present details of their results, 
in particular, the fractional contribution of this effect to the total 
absorption. 

We take into account the secondary e- pairs produced in 7- 
ray absorption by adding them to the electron injection rate. How- 
ever, our calculations underestimate the contribution of the sec- 
ondary pairs because we consider only the absorption of the ra- 
dially escaping 7-rays, for which the optical depth with respect to 
the pair production is the lowest. This approach is dictated to the 
one-dimensional nature of our numerical calculations. To take this 
process more accurately into account, three-dimensional calcula- 
tions would be required, which is beyond the scope of this work. 



3.4 Energy losses of ~10-MeV electrons 

In the spectral model above, we have identified the major radiative 
processes taking place in the source. Here, we provide estimates 
of their relative importance in various energy, density, and spatial 
regimes. We concentrate on the MeV electron energies, important 
in the radio and X-ray regimes, and on the model with the Comp- 
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Figure 4. The surface brightness profiles in three energy bands, with (a) and 
(b) corresponding to the radial profiles and spectra in Figs.[2jb),[3la) and 
Figs [2lc),[3lb), respectively. There is no 1-GHz curve in (b) because that 
model yields virtually no radio emission. 



ton emission dominating the 1 keV-10 GeV emission, and the syn- 
chrotron process producing the observed radio emission (Fig. [3^). 

Owing to the very high luminosity of the Be star, a substantial 
part of the X-ray emission can be from Compton scattering of the 
stellar radiation. At photon energies of e ~(1-10) keV, the electrons 
have 

vl/2 n „ ir,4 i^\l/2 



10 



(^— ) (j^) ^'^^^ '''' 



Electrons of about the same energy, 



(27tv) 



l/2„3/2j/2 



£1/2^1/2 



can producte GHz radio photons at S ~ 1 G. 

The high-energy electrons can be retained in the stellar-wind 
inhomogeneities, see Section 13.21 by a strong enough magnetic 
field. Assuming that the electrons diffuse in disordered magnetic 
field, we find they are retained by the magnetic field inhomo- 
geneities (in the Bohm diffusion regime) for the time. 



fdiff ■ 



2Ec 
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\10"cm/ UG/ 



B \/10MeV 
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(26) 



which is much longer than the wind t^^,^, equation illi - 

Comparing f^ifj and ^sc with the Thomson-regime IC energy 
loss time of 10 MeV electrons. 
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we see that electrons captured by the stellar wind inhomogeneities 
lose a substantial fraction of their energy before they escape from 
within the orbit. The synchrotron energy loss time is given by, 

dmnlc^ 7/lG\VlOMeV\ 
= ^ - 4 X 10' — — s, (28) 
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where equation l l25b has been used in equation i29\ . The Thomson- 
regime IC losses dominate over the synchrotron losses up to the 
(ii-independent) distance of 
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(30) 



However, this estimate is given for a fixed 6. If i? = Bq(Dq/D), 
as characteristic to stellar winds, and which was assumed in Sec- 
tion |3.3l the IC losses in the Thomson regime dominate everywhere 
(unless Bq is of the order of tens of G). 

At the density of the stellar wind, see equations ([TJ, l|4]l, 
the electrons may suffer from strong Coulomb energy losses. The 
Coulomb loss time is given by. 
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( lGouldll975l) . where /le = /liyUj/yUe. The Coulomb loss time becomes 
shorter than fic, equation J27t . below 
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(32) 



Most (but not all) of the power injected in the form of electrons with 
E < Eq goes into heating of the stellar wind (see Section|2j2j rather 
than into the synchrotron and IC emission. Equating the Coulomb 
break energy to the energy of electrons emitting X-ray IC radiation 
at e, equation ( 125 l l, we find the critical density. 



"ic=c - 4 X 10 



(3kev)(3 X lO'^cm) 



(33) 



At higher densities, the X-ray spectrum is hardened by the ef- 
fect of Coulomb losses on the electron distribution. This effect is 
only moderate in the isotropic wind, see equation Q. However, 
Coulomb losses will strongly dominate in the dense disc, see equa- 
tion O. As both ;iic=c and «„ scale as D"*, the ratio between the 
IC and Coulomb loss times is independent of the distance. 
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Then, the density above which Coulomb losses dominate over 
synchrotron losses can be found from equations ( I25t , l |29l ) and l l31b . 
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For the fast isotropic wind, equation l|4j, this density is achieved at 
a distance. 



Ds=c ^ 7 x 10 
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Again, this estimate is given for a fixed 6. If S = Bo(Do/D), the 
Coulomb losses dominate over synchrotron ones everywhere (for 
radio-emitting electrons and unless Bq is of the order of tens of G). 
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4 STRUCTURE OF THE COMPACT PULSAR WIND 
NEBULA 

VLBA monitoring over an entire orbital cycle (DMR06) have re- 
vealed an extended radio source of a variable morphology with the 
overall size of -5-10 maflat 8.3 GHz (e 3 x 10"^ eV], which 
can be identified with the compact PWN of LS I +61°303, similar 
to the compact PWN of PSR B 1259-63 jNeronov & Chernvakov3 
l2007h . At d = 2 kpc, the source size corresponds to ZJ — (1.5— 
3) X 10'"* cm. However, note that the dynamic range of the images 
in DMR06 is up to a factor of 64, i.e., the flux at the boundary of 
a radio image may be only a ~2 per cent of that at its centre. Thus, 
an e-folding distance is ~4 times less than the above overall size, 
i.e., OpwN ~ (4-8) X 10'' cm. This is in good agreement with our 
theoretical modelling of the radio flux profiles in Section [331 see 
Fig.Ha). 

For v„ ~ 10** cm s"', the wind escape time scale, Dpwn/^w ~ 

5 X 10' s. Thus, it is substantially shorter than the orbital period, 
P - 2.3 X 10* s. The strong variability of the radio morphology of 
the source over the course of an orbit can be then associated with 
this time scale. 

An interesting issue is here what is the cause of the gradient 
of the radio source structure on the scale of Opwn- It is possible 
that the decline of the flux with distance is due to the energy loss 
of electrons in the magnetic field there. Equating the synchrotron 
energy loss time, equation l |29t , to £)pwn/Vw, we obtain a relatively 
strong magnetic field, B ~ 10 G. On the other hand, this assumption 
does not need to hold, the magnetic field can be much weaker, and 
the e-folding scale may correspond to a gradient of the magnetic 
field, as we have assumed in Section [373l Consequently, the radio- 
emitting electrons lose most of the energy within the radio source 
and the radio radiative efficiency is high in the former case, or they 
lose only a small fraction of their energy and the efficiency is low 
in the latter case. 

Another important piece of information is the position of the 
peak or centroid of the radio emission (DMR06). At 8.3 GHz, the 
distance of the emission peak from the Be star varies with the or- 
bital phase in the range of £) (0.5-8) x 10'' cm, and it is always 
outside the binary orbit. At 2.3 GHz, the centroid distance equals 
D ^ (1-10) X 10" cm. This roughly agrees with our first model of 
Section [331 There is also a significant lag of the 2.3 GHz emission 
with respect to the 8.3 GHz one. Similarly to the case of the flux 
spatial gradient, the lag can be due to the electron energy loss if the 
magnetic field is high, or due to the gradient of the B field if it is 
low. 

The avoidance of the orbital region by the radio emission also 
approximately agrees with our estimates of equations U3]414b of 
that region being optically thick to free-free absorption by the stel- 
lar wind. An important unknown factor, however, is the degree of 
the heating of the stellar wind by the X-rays and y-rays, as dis- 
cussed in Section |T2l 

The alternative model of accretion/jet for the evolving ra- 
dio structure has recently been discussed in some detail by 
iMassi & Kaufman Bemaddl | |2009|) . It can explain some aspects of 
the evolution of the orbital flux periodicity, but the rotation of the 
structure with the orbital phase remains unexplained. 



On fig. 3 of DMR06, 1 cm corresponds to 6 mas, V. Dhawan, personal 
communication. 



5 THE 7-RAY LIGHT CURVE 

The IC cooling time is the shortest at the boundary between 
the Thomson and Klein-Nishina regimes for Compton scattering, 
3kT,E ~ {m^c^f, when £ ^ 30 GeV and he ^ 10(Z)/3 x 10'^ cm)^ 
s, see equation The IC scattering on TeV electrons proceeds 
in the Klein-Nishina regime. In this regime, e — E, and the electr on 
energy loss time grows with energy ( iBlumenthal & Goulj|l970h . 
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(37) 



where the density of the blackbody photons has been diluted by 
the (D/R,)^ factor. This time scale is much shorter than the energy 
loss time for 10 MeV electrons, equation ( 1271 ). The highest energy 
electrons are not retained by the stellar wind clumps because their 
diffusion time is very short, equation l l26b . Thus, the electrons emit- 
ting at the highest energies are initially present only in the vicin- 
ity of the initial injection region, around the bow-shaped average 
contact surface of the pulsar and stellar winds. Thus, the highest- 
energy emission from the pulsar/stellar wind interaction region can 
be modelled, in the first approximation, by a regular bow-shaped 
contact surface. 

Within this model, the y-ray emitting plasma flows along the 
bow-shaped contact surface of the pulsar and stellar winds, see Fig. 
[T] The contact surface is asymptotically a cone with a half-opening 
angle, 0„, which depends on the ratio of the momentum fluxes of 
the pulsar and stellar winds, 



180° 



1+7? 



(38) 



where t] 



?puisar/(A?wVwC) (R07). The shocked pulsar wind 



can flow along the cone with a large bulk Lorentz factor 
teogovalov et al.l2008h . leading to a beaming of the emission along 
the surface. Thus, a maximum of the y-ray emission is expected 
when a part of the cone passes through the line of sight. The cone 
can pass through the line of sight twice per orbit, so, in principle, 
two maxima of the y-ray light curve are expected. The true anomaly 
corresponding to those two events is at equal a ngles from the infe- 
rior conjunctioi i, cD = ± AO). T he MAGIC jAlbert et al1l2006l) 
and VERITAS jAcciari et al.l2008l) observations of the source indi- 
cate that the maximum of the TeV light curve is at the phase interval 
of 0.6 < (p < 0.7, see the solid arrow in Fig. [5] This interval corre- 
sponds to 153° < O < 167°. Taking <l>i„f = 33° (see Section [riTl 
and choosing the plus sign above, we obtain 120° < AO < 134°. 

The shift, AO, is related to Oo, and the inclination, i, via 
jNeronov & Chemvakovall2008h 



cos AO : 



cos 9„ 



(39) 



We can see that |cos^?„| < sin; is required for the cone to pass 
through the line of sight. For ; = 59° (Section [2. It and the above 
range of AO, we obtain 115° < (9oo < 127°. Interpreting this mea- 
surement of the phase of the y-ray maximum as being due to the 
passage of the cone through the line of sight, we find the ratio of 
momentum flux in the pulsar wind to that in the stellar wind of 
1.8 < 7/ < 2.4. The suggested phases of the two maxima are shown 
by the arrows in Fig.|5] 

The above simple geometrical model predicts a second maxi- 
mum of the TeV light curve, at the phase of the second passage of 
the emission cone through the line of sight. However, similarly to 
the case of the X-ray light curve, the second maximum (shown by 
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Figure 5. Attenuation of 1-TeV y-rays produced in interaction of tlie pulsar 
wind with tlie disc of Be star as a function of the orbital phase. The red 
sohd curve is for photons emitted isotropically at the position of the pulsai'. 
The blue dotted curve is for isotropic emission at D = 4R» (close to the 
boundary of the truncated disc around the Be star along the line connecting 
the two stars, see the blue box in Fig.[T). The green dashed curve is for a 
point at D = 4R» but shifted by S^b = -90° , see the green box in Fig. [T] 
The solid arrow shows the approximate observed phase of the maximum of 
the VHE y-ray light curve. The dashed arrow shows the phase of the second 
maximum, which is not observed because of the absoiption. 

the dashed arrow in Fig.|5]l is suppressed because of the compact- 
ness of the system. In the case of the TeV light curve, the suppres- 
sion is due to absorption of the emitted y-ray photons in result of 
the pair production in the photon field of the Be star. 

The effect of the yy pair production on the proper ties of the 
TeV e mission from the system was previously studied bv lSednarekl 
( l2006 a b) for the case of emission from a jet (assuming that the 
compact object is a black hole which emits jet in the direction nor- 
mal to the orbital plane) and bv lDubusI l l2006ah . under the assump- 
tion that the TeV flux comes directly from the compact object. 

In the case of the interacting pulsar and stellar winds, the re- 
gion of y-ray production is extended. Even in the simple model 
with a smooth stellar wind, see Section [JTl the y-ray emission re- 
gion extends along the bow-shaped interaction surface. Thus, ab- 
sorption of the y-ray emission emitted toward an observer varies 
strongly across the surface. The strongest absorption is for photons 
emitted around the apex of the bow-shaped surface, closest to the 
massive star. The situation becomes even more complicated when 
the stellar wind is clumpy. Thus, we can estimate only limits within 
which the optical depth, Tyy, varies across the emission region at 
different orbital phases. 

The results of such calculations are shown in Fig. [5] which 
takes into account the anisotropy of UV emission from Be star. The 
red solid curve shows the attenuation of 1-TeV photons emitted at 
the position of the compact object. This point can be considered as a 
rough outer boundary of the pulsar/stellar wind interaction region. 
(Note t hat the shape o f the curve slightly differs from that calcu- 
lated bv lDubusll20063. because we use the current estimates of the 
orbital parameters of G07, whereas he used the results of COS.) The 
blue dotted curve shows the attenuation of photons emitted at the 
point situated at the distance of D = 4^, 2 x 10'^ cm from the 
centre of the Be star on the line connecting it to the compact ob- 
ject. This distance is a rough estimate of the radius of the truncated 
Be-star disc, see Section [2^ and we also used it in our spectral 



calculations in Section [331 Finally, the green dashed curve shows 
the attenuation of emission produced at Z) = 4if,, but <5<I> = -90° 
away from the direction of the line connecting the two stars. These 
three points at = 0.7 are shown by the red, blue and green box, 
respectively, in Fig.[T] 

We see that the orbital dependence of the attenuation varies 
by a large factor across the y-ray emission region. The transmitted 
fraction reaches the maximum at the inferior conjunction for the 
emission on the line connecting the two stars, but that maximum 
is shifted to a later phase for the blue dotted curve. An accurate 
calculation needs to assume the detailed dependence of the emis- 
sivity on both the position and the phase. This cannot be done as 
the morphology of such extended sources remains so far unknown. 



6 CONCLUSIONS 

We have shown that the compact PWN model developed here can 
explain a number of observational properties of LS I -l-61°303, from 
the radio to the TeV y-ray energy band. The main new feature of 
this model is to account for the dumpiness of the fast wind com- 
ponent from the Be star, which slows down the escape of the high- 
energy electrons from the system. Also, due to plasma instabili- 
ties, the relativistic electrons are retained in the moving clumps 
of the stellar wind and carried away with them. The presence of 
time-dependent X-ray emitting clumps explains the X-ray variabil- 
ity observed on time scales much shorter than the orbital period. 
On the other hand, our model has not resolved yet the issue of the 
momentum flux of the pulsar wind being significantly higher than 
that of the Be wind, which presents a problem for interpretation of 
the observed radio structures (as pointed out by ROT). 

The electrons lose energy in the IC, synchrotron and Coulomb 
processes. We have presented detailed inhomogeneous spectral 
models, reproducing the average spectra of the object. The 1 keV- 
10 GeV spectral range is dominated either by synchrotron or IC 
emission, depending on the form of the electron injection. Our 
models reproduce very well, in particular, the sharp spectral cut- 
off seen in the GeV range by Fermi. 

We have also provided detailed formulae and calculations 
comparing the relative importance of those processes as functions 
of the electron energy and the position in the source, measured by 
the distance from the Be star. Coulomb energy losses are impor- 
tant for MeV electrons. The radio emission is suppressed by strong 
free-free absorption, up to the distance of about an AU from the 
source, in agreement with the radio data. The free-free absorption 
also leads to a strong attenuation of the pulsed radio emission from 
the neutron star. Compton losses dominate in an inner part of the 
fast wind. We also find that a part of the wind surrounding the high- 
energy source is strongly heated by both Compton and Coulomb 
processes. 

We find the most likely mechanism explaining the observed 
orbital modulation at TeV energies is anisotropy of emission, with 
relativistic electrons accelerated along the surface at which the 
pressure of the pulsar wind equals the average pressure of the stel- 
lar wind. In this model, pair absorption of the TeV emission plays 
a role by suppressing the one of the two expected maxima. 
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